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ABSTRACT 

Structure and spectral energy distributions (SEDs) of externally irradiated circumstellar disks are 
often computed on the basis of the two-temperature model of Chiang & Goldreich. We refine these 
calculations by using a more realistic temperature profile which is continuous at all optical depths and 
thus goes beyond the two-temperature model. It is based on the approximate solution of the radia- 
tion transfer in the disk obtained from the frequency-integrated moment equations in the Eddington 
approximation. We come up with a simple procedure ( "constant g z approximation" ) for treating 
the vertical structure of the disk in regions where its optical depth to stellar radiation is high. This 
allows us to obtain expressions for the vertical profiles of density and pressure at every point in the 
disk and to determine the shape of its surface. Armed with these analytical results we calculate the 
full radial structure of the disk and demonstrate that it favorably agrees with the results of direct 
numerical calculations. We also describe a simple and efficient way of the SED calculation based 
on our adopted temperature profile. Resulting spectra provide very good match (especially at short 
wavelengths) to the results of more detailed (but also more time-consuming) SED calculations solving 
the full frequency- and angle-dependent radiation transfer within the disk. 

Subject headings: accretion, accretion disks — circumstellar matter — stars: pre-main-sequence 



1. introduction 

Protoplanetary disks around T Tauri stars (stellar 
mass M* < 2 M©) as well as the disks around more 
massive Herbig Ae/Be stars (M* > 2 M©) belong to the 
class of circumstellar disks that are strongly affected by 
the radiation of their central stars. Irradiation modifies 
both the vertical structure of the disk and the radial de- 
pendence of disk properties. As a result, observational 
signatures of centrally irradiated discs are quite different 
from those of the conventional viscously heated a-disks 
(Shakura & Sunyaev 1972). 

Kenyon & Hartmann (1987) have pointed out that the 
observed infrared spectra of protoplanetary disks can be 
understood as resulting from passive reprocessing of the 
stellar radiation assuming that disk surfaces have flared 
geometry. Later Chiang & Goldreich (1997; hereafter 
CG97) have analytically derived hydrostatic, radiative 
equilibrium models of disks around T Tauri stars. Their 
calculation was based on a two-temperature approxima- 
tion for the disk thermal structure in which the outer 
"superheated" layer has been assumed isothermal at the 
temperature characteristic for dust grains directly illumi- 
nated by unattenuated stellar radiation. The isothermal 
midplane region was found to have lower temperature set 
by the balance between the oblique stellar illumination of 
the disk surface and the isotropic infrared emission of the 
disk interior. The success of this simple model has led 
to its wide acceptance for interpreting the observations 
of circumstellar disks. 
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The goal of this paper is to improve our understanding 
of the circumstellar disk structure by going beyond the 
simple two-temperature approximation of CG97. Cal- 
vet et al. (1991; hereafter C91) have studied the vertical 
radiative transfer within irradiated dusty circumstellar 
disks and obtained approximate analytical expression for 
the vertical temperature profile which is superior to the 
two-temperature approximation of CG97. In this work 
we employ this temperature distribution, which is re- 
viewed in detail in <j3J to calculate the disk structure 
and spectrum. Under certain conditions it is possible to 
derive approximate analytical solutions for the vertical 
disk structure as demonstrated in fJ3 This allows us to 
efficiently compute the radial dependence of various disk 
properties in <21 an d to compare them with purely nu- 
merical solutions. Finally, in SjSlwe demonstrate how the 
use of the realistic temperature profile of C91 improves 
the calculation of the spectral energy distribution (SED) 
of circumstellar disks. 

Throughout this study we assume the major heating 
source of the disk to be the radiation of the central star 
of mass M*, having radius R* and effective temperature 
Viscous heating is considered subdominant. Circum- 
stellar disk extending from Ri„ to R ou t is assumed to 
be geometrically thin thus allowing us to neglect radial 
transport of radiation. This simplifies radiative trans- 
fer within the disk and reduces it to a one-dimensional 
problem. Surface density of the disk So varies as 

£ (a) = £* , (1) 

where a is the distance from the central star and X* = 
E(i?*). We use 5 = 1 in our calculations. Our re- 
sults will be illustrated using three fiducial models of 
circumstellar disks around Herbig Ae/Be star, T Tauri 
star, and a brown dwarf, analogous to those adopted by 
Dullemond & Natta (2003; hereafter DN03). Proper- 
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ties of these models are summarized in Table IHT1 (£] = 
AU) 5 is the value of S at 1 AU). 

2. THERMAL STRUCTURE. 

To describe the transfer of stellar radiation with black- 
body temperature within the disk we introduce optical 
depth t*(z) (z is the altitude above the disk midplanc) 
which is calculated along the radial direction from the 
star. Following conventional wisdom we approximate t* 



10 2 10 3 



10 5 T, K 



K P (T*)p(z)dz = a- 1 K P (T*)S(z), (2) 



where p(z) is the local gas density, S(z) = p(z)dz is 
the disk surface density above a given z, np(T) is the 
Planck mean opacity at temperature T, and a is the 
flaring angle (tt/2 — a is the angle between the normal 
to the disk surface and the radial direction). In the thin 
disk approximation 
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where Hi (a) is the altitude of the disk surface defined 
as a surface at which T*(a, Hi) = 1 and r\ « 0.4 is con- 
stant (Kenyon & Hartmann 1987; CG97). There are two 
regimes of irradiation: "lamp-post" illumination when 
R* ^ Hi and 



R* 
v — , 

a 



and "point-source" illumination when R+ < Hi and 

d (Hi 
a « a— ' 

da 



(4) 
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In the first case the geometry of illumination is deter- 
mined by the finite size of the source, while in the sec- 
ond it is the self-adjustment of the disk surface geometry 
that sets the illumination angle. Apparently, disk must 
be flared in the latter case, i.e. d{H\/a)/da > 0. 

Dust grains are the major opacity source in circum- 
stellar disks except their central regions where high tem- 
perature leads to dust sublimation. To compute Kp(T) 
we use the frequency-dependent opacity k v characteris- 
tic of 0.1 /im silicate grains (Draine & Lee 1984) with no 
scattering (the same as k v used by DN03 6 , to facilitate 
subsequent comparison with their results). Both Kp(T) 
and k v used in this study are shown in Figure ^ 

At low temperatures (T < 10 2 K) up behaves as a 
power law of T. Motivated by this, in deriving our an- 
alytical results we will frequently use the "local power 
law opacity" approximation to Kp(T). It is defined as 
the best power law fit to the actual run of Kp(T) in the 
temperature range relevant at a given distance from the 
star a: 

n P {T) a K (a)T^ a \ (6) 

where parameters Ko and (3 are functions of a. This sim- 
plified approach shows very good agreement with numer- 
ical calculations using full Kp(T), see 21 

6 See http://www.mpia-hd.mpg.de/homes/dullemon/radtrans/ 
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Fig. 1. — Opacity k v (solid curve) adopted in this paper as a 
function of frequency v (lower axis) and Planck opacity Kp (dashed 
curve) computed based on this function of temperature 

T = hv/k (upper axis). 



Calvet et al. (C91) have derived approximate tem- 
perature structure of centrally irradiated accretion disk 
heated from inside by viscous dissipation. Their result is 
based on the assumption that disk is optically thick to its 
own radiation with opacity independent of temperature. 
For simplicity, in this study we neglect radiation scat- 
tering by dust and viscous heating. At the same time, 
the assumption of a geometrically thin disk allows us to 
generalize C91 result for an arbitrary dependence of Kp 
on T . Simple calculation in the spirit of C91 results in 
the following implicit expression for T: 



T\z) = ^ 



a 



a4>ex g(T) (z) 
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(7) 



Here q(T) = Kp(T i )/ Kp(T) is a ratio of opacities at tem- 
peratures and T, factor (f> represents the fraction of 
the stellar surface visible from a given point of the disk 
surface, while r*(z) and a are given by equations (J2J & 
©. If the inner edge of the disk is close to the stellar 
surface - an assumption that we adopt in this paper - 
then 4> sa 1/2. Following Dullemond, Dominik & Natta 
(2001; hereafter DDN), we have introduced correction 
factors ipex < 1 aud ?psh < 1 to make this solution valid 
even when the disk is optically thin to the radiation of 
its outer layer and/or its interior [situations not allowed 
in the original optically thick solution of C91 which can 
be recovered from eq. Q by setting ip e x = ipsh = 1]- 
We provide details of the calculation of -0-factors in Ap- 
pendix 

According to equation JJJ the vertical extent of the disk 
can be split into two layers: the outer layer exposed to 
stellar radiation and the inner layer shielded from direct 
starlight (see Figure for an illustration). 

Dust within the exposed layer absorbs highly 
anisotropic stellar flux with blackbody temperature T* 



3 



Exposed 
layer 




Ho Hi 



Fig. 2. — Schematic illustration of the vertical thermal struc- 
ture of irradiated disk adopted in this paper. See text for detailed 
explanations. 



and reradiates it isotropically at lower temperature. A 
half of the reradiated flux escapes while another half il- 
luminates the interior of the disk. Temperature in the 
exposed layer varies vertically according to 

1/2 
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where 



T ex (a) = T4MTe X )} 1/4 [ ^ 
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Expressions © and (jl 1|> correspond to the case of the 
power law opacity I©- It follows from equation JHJ that 
T ex is the temperature of the optically thin part of the 
disk (t* <C 1). Exposed layer corresponds to < t c 
where 



= In 



q{T S h) ipsh 



(12) 
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We denote H (a) = z(t+ 
cal depth t* equals r c . 

Shielded layer lies below the exposed region at high op- 
tical depth, t* > r c . Within this layer disk material re- 
ceives only the reprocessed radiation of the exposed layer 



since the direct stellar flux is almost entirely absorbed 
by the exposed layer. Shielded region is isothermal with 
temperature 



T sh (a) = T* 



2 ipsh 



1/4 



a 



1/2 



n > t c . (13) 



Equation (|13Jl results from balancing the incoming flux 
of the reprocessed radiation with the energy loss from the 
shielded layer. Factors ip s h and ip ex in equations l(12|) and 
1)13(1 accounting for the possibility of optic ally t hin disk 
are set by T s h and T ex , see equations IJA2I) and (|A4|> . 

Temperature profile J7J should be compared with the 
temperature structure adopted by CG97: 

T = T ex , t* <C 1 , 
= T sh , r* > 1. (14) 

Equation provides a better match for the real thermal 
structure of the exposed layer as it does not feature an 
unphysical temperature jump at some optical depth. 

3. VERTICAL STRUCTURE. 

Hydrostatic equilibrium in z-direction is described by 
dP 

dz 

where P is the gas pressure and = (GAf^/a 3 ) 1 / 2 is the 
local angular frequency. Within the isothermal shielded 
layer (z < Hq) this equation yields 



-tfzp, 
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p(a,z) P(a,z) 



exp 
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po(a) P (a) 

where h s h = fl~ 1 (kBT s i l /p) 1 / 2 is the isothermal scale- 
height of the disk in the shielded region, and po and Pq 
are the midplane values of gas density and pressure. Be- 
cause of the exponential decay of p with z shielded region 
contains most of the surface density So (a) of the disk (if 
Hq > h s h) which implies that 

-1/2 S o( a ) 



Po (a) ~ (27r) 
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3.1. Constant g z approximation. 



Exposed region is not isothermal when ~ 1 but we 
can still understand its vertical structure by introducing 
the so-called constant g z approximation. This approxi- 
mation makes use of the fact that gas density in the disk 
decreases very rapidly with increasing z, see e.g. (116(1 . 
As a result, when considering the hydrostatic equilib- 
rium, we can to zeroth order neglect the variation of the 
vertical acceleration g z = fl 2 z compared to the change 
of p with z in equation l|15l) . 

In particular, for the determination of the disk struc- 
ture near its surface we can set z ~ Hi which allows 
us to integrate equation (|15fl over z. We find that 
P(z) w Q 2 Hi'E,(z) and using definition © obtain the 
following relation between P and r*(z): 



n 2 H 1 

a tt^tT*. 



Substituting this into equation l|15f) one finds 

1 dr+ £l 2 p 
r* dz k B T(T+)' 



(18) 



(19) 
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In the shielded region below Ho (we assume that Ho ~ 
H\ which is verified below) T(t*) ~ T s h and equation 
P|) yields 



P(z) ^ p(z) 



2/i 2 



In 



- hir*. 



(20) 



Integration constant is fixed in equation 1(201) with the aid 
of equation (|18fl by noticing that in the shielded region 
solution (|2"U)l must reduce to (fTSj) . Substituting = t c 
into equation 1)20(1 one obtains the height of the boundary 
between the exposed and shielded layers Hq'. 
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Pqk(T*) 
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= (21nA) 1/2 
Sok(T + ) 



(21) 
(22) 



Our use of Hq instead of Hi in the definition of A im- 
proves the accuracy of ill' 21) . 

Structure of the exposed region will be determined 
under the assumption of the "local power law opacity" 
approximation © which allows us to use equation @. 
When integrating equation 1(19(1 for z > Hq we can split 
the integral into two parts: J * = j c +/ *• In the first 
integral we set T(t+) = T s h which reduces it to equation 
(J5DJ), while in the second we use T(t*) given by equation 
(0. As a result, for z > H one finds 

22 ^»- 2 ^[ Ei (-lf?)- Ei (-IT?)]< 23) 

where h ex = fl~ 1 (ltBT ex / p,) 1 / 2 is the disk scale 
height corresponding to temperature T ex and Ei(x) = 
— t~ 1 e~ t dt is exponential integral (Gradshtein & 

Ryzhik 2000). 

Equation ((23(1 allows us to determine the position of 
the disk surface: 



Hl^Hl + X 2 h 2 exl 
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(24) 
> 0, (25) 



where for all practical purposes factor \ ~ 1 can be con- 
sidered almost constant - although r c > 1 is a function 
of a, the corresponding contribution in i|25|) is small com- 
pared to the leading term Ei(— 1/(4 + (3)). One can see 
that Hi r*> Hq unless A < xT ex /T s h (which can only 
be realized when the constant g z approximation breaks 
down, see H4 4. II and Figures OEJ ■ Asymptotically 

Ei{-x)~-—\l + 0(x- 1 )} 1 x>l, 
x 

~C + ln(x) +0(x), x<l, (26) 

(Gradshtein & Ryzhik 2000) where C » 0.5772 is the 
Euler constant. These limiting forms suggest that 

-n/(4+/3) 



z 2 ^H 2 + 2(4 + fi)h 
forl<r* <r c 7 (£T < z < H x ) 
t*(z) w (4 + /3)e- c exp 



2 ' 

ex 



while 

z 2 - 
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(28) 



7 This assumes that r c ^> 1 which may be a questionable as- 
sumption in real disks. 
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for < 1 (z > i/i). The dependence (|28|l is some- 
what unexpected given that the exposed region above Hi 
is virtually isothermal at temperature T ex . Apparently, 
this is an artifact of our adopted constant g z approxima- 
tion. At the same time P(z) and p(z) given by l|29|l are 
in accord with the isothermality of this part of the disk 
which is to some extent coincidental. 

Equations (|21|I - H29(I fully determine the vertical hydro- 
static structure of the disk as well as the positions of 
Hq (a) and Hi (a) through the parameters of the disk and 
the star. 

3.2. Comparison with the two-temperature 
approximation. 

It is instructive to compare the results of the previ- 
ous section with the disk structure according to CG97. 
Their model exhibits discontinuous jump of temperature 
at t* = 1 and assumes upper disk layer (t* < 1) to 
be fully isothermal, see equation 1(14(1 . This layer ("su- 
perheated dust layer" in their nomenclature) plays the 
role of our exposed region, while the underlying region 
(t* > 1) should be identified with our shielded layer. In 
both layers vertical profile of pressure is Gaussian and p 
is discontinuous across the boundary defined by t* = 1 
(while P is continuous there). Apparently, H = Hi 
in the CG97 model; applying constant g z approximation 
with their temperature structure 1(14(1 one finds Hi to be 
given by equation ((21(1 but with A replaced by 



A 



CG 



Sqk(T^) 
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(30) 
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The position of the disk surface computed in the two- 
temperature approximation deviates from that given by 
equation ((24(1 . Difference is not very significant in the 
inner disk where X 2 h 2 h 3> X 2 h1x- However, at intermedi- 
ate 8 radii this is no longer valid and the two-temperature 
model underestimates Hi compared to our more accurate 
expression (J23J. 

4. RADIAL PROPERTIES OF IRRADIATED DISKS. 

Based on analytical results obtained in the previous 
section we can provide a complete description of the ir- 
radiated disk structure. This requires the determination 
of T s h(a) which can be done by substituting -Hi (a) given 
by the expression ((24(1 into © , plugging a into ((13(1 and 
solving the resultant differential equation for T s h(a). 

Analytical expressions for T s h(a) can be obtained for 
optically thick disks (ip ex = ip s h = 1) which either (1) 
are illuminated in the lamp-post regime or (2) experi- 
ence point-like illumination and have A = Ho/h s h of the 
order of several (typically > 3). In the latter case one 
can neglect in equation ((24(1 term xhex compared to Xh s h, 
which, coupled with the weak sensitivity of A to a, leads 
to Hi oc h s h- This assumption should be valid in the in- 
ner, dense parts of circumstellar disks. Analytical results 

8 At large r adii ne ither CG97 nor constant g z approximations 
work well, see i|4 4.11 
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Fig. 3. — Comparison between the semi-analytical (constant g z 
approximation) and numerical approaches to computing the disk 
structure in Model I. (a) Position of the disk surface H\ (a) /a 
calculated using eq. 1241 (solid curve) and numerically with 
the full treatment of the vertical hydrostatic equilibrium (short- 
dashed curve). Also shown are the semi-analytical results for 
Hg(a)/a (dot-dashed line), h a f l (a)/a (long-dashed line), and the 
ratio hsh/Ho = A -1 (dotted line), (b) Midplane temperature 
Tsh( a ) computed semi-analytically (solid curve) and numerically 
(short- dashed curve). Two methods give virtually identical results 
in this case. Dotted line represents T ex (a). (c) Flaring angle ct(a) 
obtained semi-analytically (solid curve) and numerically (short- 
dashed curve). 



for these two cases are summarized in Appendix [H] In 
particular, one finds that Hi oc a 9 / 8 for lamp-post illumi- 
nation and Hi oc a 9 / 7 for point-like illumination, i.e. in 
both cases disk surface is flared. In a more general case 
we find the behavior of T s h(a) by solving the system of 
equations 10, l|13[l. and lj2"4T) numerically. 

To test the validity of the constant g z approximation 
we additionally compute disk structure without analyt- 
ical approximations of i J3 3.11 with the only assumption 
that the temperature profile is given by equation J7J. In 
this approach at every radius we calculate vertical disk 
structure numerically with eq uation s Q) and (|T3|l . in- 
stead of using the results of W3 3.11 The value of Hi 
obtained in this fashion (instead of using equation |24p 
is then used in equation (J3J) to reconstruct the radial disk 
profile. This approach uses the detailed opacity descrip- 
tion Kp(T) shown in Figure ^ while our semi- analytical 
procedure calculates disk structure assuming local power 
law opacity fit © to the curve of Kp(T). Some details 
of the numerical disk structure calculation can be found 
in Appendix IO 

In Figures T3I5I we pre sent our results for the three disk 
models listed in Table O We compare the behaviors 
of -Hi (a), T s h(a), and a(a) found semi-analytically and 
numerically. One can see that the theory based on the 
constant g z approximation predicts -Hi (a) quite well - 
the discrepancy with the numerical result at the level 
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Fig. 4. — Same as Figure|3]but for Model II. 



of ~ 10% appears only in the outermost parts of the 
disk. The same is true for T s ^(a). Flaring angle a is a 
more sensitive probe of the difference between the two 
approaches as it involves a derivative of Hi(a). In the 
outer regions of disks in models II and III semi-analytical 
a(a) diverges significantly from the numerical solution. 
This is the result of the lower optical depth in these mod- 
els compared to the model I, which affects the validity of 
the constant g z approximation, see H44.ll 

Nevertheless, at small and intermediate a the overall 
agreement between the two methods of the disk struc- 
ture calculation is very good. For illustration, in Figure 
E|we present the comparison between the vertical profiles 
of T, P, p, and r obtained by both methods at a = 10 
AU in Model II. Despite the fact that at this location 
constant g z approximation is already close to the limit 
of its validity (see Figure QJ the agreement in the ver- 
tical structure of P, p, and T is quite remarkable. The 
constant g z approximation fails only for t*(z) which was 
expected, see discussion after equation (|29Jl . 

4.1. Validity of the constant g z approximation. 

Constant g z approximation should work well when 
Ho, Hi > h s h since in this case the decay of P and p 
with height is well represented by the tail of the Gaus- 
sian profile at z ~ Ho, Hi. As a result, at this height p 
varies with z much faster than g z and constant g z limit 
works well. Based on this reasoning, we can tentatively 
suggest Hq > 2h s h as an approximate criterion for the 
validity of this approximation. 

This conclusion is confirmed by examination of Fig- 
ures 13151 which demonstrate that the relative difference 
between the values of a computed with and without con- 
stant g z approximation starts to exceed 10% when the ra- 
tio Ho/h s h = A drops below 2.5. Using definition l|21|l 
we can then claim this approximation to be accurate (at 
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Fig. 5. — Same as Figure |3] but for Model III. 

the level of ~ 10%) whenever A > 20 or 

KP (T.)E > 10, (31) 

where we have taken into account that in our models 
\/27rar c ss 0.5 at the limit of applicability. This condi- 
tion is accurate only in the order of magnitude as the 
exact factor in its right hand side is exponentially sensi- 
tive to the uncertainties in A. Despite this, in practical 
situations one may still use condition 1)3 l|l as a rough 
proxy for checking whether the constant g z approxima- 
tion can be relied upon. 

5. SPECTRAL ENERGY DISTRIBUTION. 

In this section we calculate disk SED using our adopted 
temperature profile (JJJ, and then compare the outcome 
with the results of other approaches. 

In our case SED is produced by emission of both 
shielded and exposed layers. The contribution of the 
former per unit surface area of the disk dF^ h {a)/dS is 
given simply by (see DDN) 

dF* h {a) ., 

, — . , , ex p 



dS 



■ 2tt cos i 
x B, 



1 



Y, (a)K v 



ATsh(a)}, (32) 

where i is the disk inclination (we assume disk thickness 
to be so small that all elements of its surface can be 
characterized by a single value of inclination), B„(T) is a 
Planck function, and the factor in parentheses accounts 
for the possibility of the disk being optically thin at a 
given frequency v. 

Exposed layer is always optically thin to its own radi- 
ation and its contribution to the SED is given by 



dS 



2tt < 1 + exp 



E (o)k v 



COS I 



B v [T(a, t*)] dr v 



(33) 
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Fig. 6. — Vertical structure of the disk in Model II at a = 
10 AU. Results for p(z) (a), t(z) (5), T{z) (c), and P(z) (d) are 
shown. Solid curves represent fully numerical results while dashed 
curves represent results obtained using constant g z approximation. 
Locations of h^, Ho, and Hi are marked on each panel by dashed 
lines. 
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Fig. 7. — Spectral energy distribution of a disk annulus around 
Herbig Ae/Be star (T* = 10 4 K, ij* = 2 Rq, M* = 2 Mq). 
Annulus is placed at a = 1 AU, it has a width Aa = 0.01 AU, 
flaring angle a = 0.05, and vertical optical depth r = 10 at 550 
nm. Das hed line shows SED calculated according to equations 1321 
and 1341 . solid line represents SED calculated using full frequency- 
and angle-dependent radiation transfer (DN03), while dotted line 
is SED computed in the framework of the two-temperature approx- 
imation (CG97). 
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A, fj.m A, yU.m 



Fig. 8. — Same as Figure 171 but for annulus at a = 50 AU from 
the star having width Aa = 0.5 AU and vertical optical depth 
t = 300 at 550 nm. 



Fig. 9. — Same as Figure|jjbut for an annulus around T Tauri 
star (T* = 4000 K, R+ = 2 R e , M* = 0.5 Rq). Annulus is placed 
at a = 0.1 AU and has width Aa = 0.001 AU and vertical optical 
depth t — 10 at 550 nm. 



where t„ = ar^K^/ Kp(T+) is the optical depth at fre- 
quency v along the disk normal, T(a, T*) is given by 
and the factor in parentheses accounts for the contribu- 
tion of the exposed layer on the other side of the disk 
when the shielded layer is optically thin at frequency v. 
Here r„. c = oit c k v / 'kp(T^) <C 1 is the value of r v at the 
boundary between the shielded and exposed layers. 

Using equation (|10f> one can switch from integration 
over dr v to integration over y = hv/ksT in the range 
Ho < y < oo, where yo = hv/k,BT ex . The upper limit 
hv/ksTsh can be extended to oo since most of the ex- 
posed layer emissivity originates at temperatures ~ T ex 
corresponding to y ~ yo- As a result, one obtains 



dF™(a) 



dS 



2tt < 1 + exp 



cos i 



A„(T ex ),(34) 



where 



Kp(l+) J I 




dy ip (hv/ky) 



c 2 Kp(T^) J y e y - 1 



(35) 



is the emissivity of the exposed layer. Function ip(T) is 
defined as 



^(T)=4 



din Kp(T) 



d\nT 

Jq 00 B v (T)k v (4 + din K„/dln v) dv 
J °° B v {T)n v dv 



(36) 



Equations l|34|) - (|36|l provide us with the method of cal- 
culation of the SED of the exposed layer that does not 
require an explicit solution of Q to obtain T(t*) but 



at the same time fully accounts for the non-isothermal 
nature of the exposed layer. 

Full emissivity is given by the sum of dF^ h (a)/dS and 
dF^{a)/dS. In Figures 17151 we display the SED pro- 
duced by the thin disk annulus placed at different radii 
around different stars (their parameters can be found 
in figure captions), assuming a fixed value of the flar- 
ing angle a — 0.05 and inclination i = 0° for all of them. 
This choice allows us to directly check our results against 
the calculations of DN03 who have made an exhaustive 
comparison between the two-temperature model and the 
full numerical treatment of the radiative and hydrostatic 
equilibrium of the disk, using both frequency- and angle- 
dependent radiation transfer in z-direction. 

Figures T7I9I show that at short wavelengths our treat- 
ment (|34|) - (l3fill of the temperature structure of the ex- 
posed layer reproduces full numerical SED remarkably 
well, much better than the two-temperature approxima- 
tion. This is not surprising since unlike our equation Ijl0|l 
the two-temperature model does not capture the smooth 
transition between T ex and T s h in the exposed layer and 
this affects its performance at high frequencies. 

At longer wavelengths situation is more complicated. 
In particular cases represented in Figures an d El there 
is very nice agreement between our SED calculations and 
the numerical results of DN03 in the whole range of wave- 
lengths. However, in the case displayed in Figure [S] situ- 
ation is not so simple: both our approach and the two- 
temperature approximation give very similar results but 
deviate significantly from the SED computed by DN03 
at A ~ 100 /im. This has to do with the fact that at 
longer wavelengths disk spectrum becomes more sensi- 
tive to the details of the thermal structure of the shielded 
layer which is treated similarly by both our method and 
the two-temperature approach of CG97. However, it has 
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been previously shown by Dullemond, van Zadelhoff & 
Natta (2002) that in some cases neither equation Q nor 
l|14l) reproduce the temperature structure of the shielded 
layer accurately enough, and this causes discrepancies at 
long wavelengths in spectral calculations. Nevertheless, 
the much better agreement at short wavelengths obtained 
with our procedure at the expense of only slightly in- 
creased (compared to the two-temperature model) com- 
plexity justifies its practical use in various applications. 

6. DISCUSSION AND SUMMARY. 

We have explored structure of the irradiated circum- 
stellar disk based on a realistic temperature profile Q, 
which is more accurate than the two-temperature ap- 
proximation of CG97. Our results on the disk SEDs pre- 
sented in persuasively demonstrate that this approach 
is superior to that of CG97 or DDN as it naturally repro- 
duces the spectral behavior at short wavelengths without 
making any artificial assumptions. 

This improvement is significant because of the recently 
emerged interest to the structure of the SED around 
A ~ 2 /^m. Many circumstellar disks exhibit excess emis- 
sion in this band (Hillenbrand ct al. 1992) which has been 
generally interpreted as the evidence for the existence of 
the dust sublimation region in the inner parts of these 
disks (DDN). Examination of Figures & |H1 reveals that 
the use of the two-temperature approximation can signif- 
icantly overestimate the contribution of the disk to the 
total flux in this band. As a result, one would under- 
estimate the amount of emission produced by the dust 
sublimation region and this can seriously affect the in- 
terpretation of the data. On the contrary, our approach 
reproduces realistic disk SEDs (obtained by DN03 using 
rather time-consuming procedure for solving full radia- 
tion transfer) very accurately at short wavelengths and, 
thus, can be relied upon when determining the emissiv- 
ity of the dust sublimation region from the data. Note 
that our method of SED calculation based on equation 



J7J is no more computationally intensive than the two- 
temperature approximation of CG97 which is routinely 
used for analyzing protoplanetary disk spectra. 

Another goal of this study was to introduce and test 
the constant g z approximation ( t!3 3.1 \ for treatment of 
the irradiated circumstellar disk structure. We have 
demonstrated by comparison with numerical results 
that within its region of applicability [see equation (13 II) ] 
this approximation works very well. Its use has allowed 
us to obtain analytical solutions for the vertical disk 
structure and to determine the position of the disk sur- 
face as a function of distance from the star rather accu- 
rately through the disk and stellar parameters. Improved 
description of the density and temperature structure of 
gas in the disk photosphere at t* ~ r c obtained in H3 3.11 
allows one to improve the predictions of molecular line 
intensities compared to the two-temperature approxima- 
tion [see Dullemond et al. (2002) for comparison of the 
line intensities obtained using full radiation transfer and 
its simplified representation given by equation lf7|]. 

Finally, solutions presented in A3 3. II allow fast and 
efficient exploration of the large phase space of circum- 
stellar disk parameters when fitting disk SEDs. Their 
analytical nature will also be useful for treating more 
complicated problems, e.g. determining the structure of 
the dust sublimation region, ice lines, and so on. 



RRR thankfully acknowledges the financial support 
by the Canada Research Chairs program, W. M. Keck 
Foundation, and NSF via grant PHY-0070928. FDC ac- 
knowledges financial support through the fellowship of 
the DGEP-UNAM, and from the EU Community under 
the Marie Curie Research Training Network " JETSET". 
FDC appreciates the hospitality and financial support of 
the Institute for Advanced Study (IAS) during the visit 
when this work has been initiated. 



APPENDIX 

■0-FACTORS FOR OPTICALLY THIN DISKS. 

Distant parts of the circumstellar disk where Eo is small and temperature is low can be optically thin to the radiation 
of the shielded layer (T, np(T s h) < 1) and/or the exposed layer {Y> np(T ex ) < I). To account for this possibility one 
introduces factor ip s h < 1 defined as the ratio of the shielded layer emissivity to its value in the optically thick case. 
Factor 4>ex < 1 is defined as the fraction of radiation emitted by the exposed layer that gets absorbed by the shielded 
layer. 

There are different ways of defining -0-factors. In particular, Chiang et al. (2001) used simply 
i>^h( T sh) = 1 - exp [-Y, K P (T sh )} , ip° x (T ex ) = 1 - cxp [-'S Kp(T ex )] , 
while DDN suggested more accurate expressions for ip s h and ip ex : 



(Al) 



Vsh { 1 sh) 



J °° B v (T sh ) [1 - exp (-£()«!/)] dv 



J °° B v {T ex )K u [I - exp (-£ Ki/)] dv 



(A2) 



(A3) 



Kp(T ex ) 

All of these expressions assume that corresponding layers of the disk emit pure blackbody radiation which is a 
reasonable assumption for the isothermal shielded layer (and so ipf™ should be accurate). However, the temperature 
of the exposed layer varies quite dramatically when I < < r c and its spectrum must deviate from B v (T ex ). For this 
reason, for ip ex we use the following expression: 

/ °° A v {Tex) [1 - exp (— SoKv)] dv 



J™A v (T ex )dv 



(A4) 
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Fig. A10. — (a,c) Plots of i/j ex (solid curve), ipg x (dotted curve), and V'SP N (dashed curve) for two different values of disk surface density: 
£ = 0.3 g cm~ 2 (a) anf Eq = 1.0 g cm -2 (c). (b,d) The same for j/)JL (dotted curve), and i/>J]P N (dashed curve). 



where A V (T) is the true emissivity of the exposed layer given by equation (|35|l . Clearly all i/'-factors reduce to unity 
in the optically thick disk (So k p(7 , s /i), ^oKp(T ex ) > 1). 

In Figu re |A 101 we plot different ^-factors for two values of £. One can see th at al l three prescriptions ilAlf) . I|A3[) . 
and i|A4(l for ip ex give rather similar results. At the sa me t ime, simple expression (|A1|) for ipg^ advocated by Chiang et 
al. (2001) differs from the more accurate prescription J^^J) at t he le vel of tens of per cent in some temperature ranges. 
In this work we use ipsh and ip ex given by equations (|A2|) and l|A4fl correspondingly. 
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Table CI. Representative models of circumstellar disks. 



Model 


M* [Mq] 


R* [R@] 


1* [K] 


Rin [AU] 


Rout [AU] 


Si [g cm 2 ] 


Object 


I 


2 


2 


10000 


1 


300 


1000 


Herbig Ae/Be 


11 


0.5 


2 


4000 


0.1 


300 


100 


T Tauri 


111 


0.1 


1.3 


2600 


0.033 


30 


100 


Brown Dwarf 



RADIAL PROPERTIES OF IRRADIATED DISKS. 

In the case of lamp-post illumination one immediately finds from |j2J and (|13|l that 

^°> = r -(ff (tP «">^(?f (if • < B1 > 

where h* = Qr^fcT*/^) 1 / 2 , = (GM±/Rl) 1/2 . Value of A in this case is given by 

^"p^fels:) • (B2) 

In the case of point-like illumination the substitution of equation (J5J with Hi « _ff = Afo s /j m to definition l|13|l 
yields a differential equation for T s /,., which can be easily solved if A is roughly constant. One finds 

For point-source illumination A is given by 

ps ~ (87r) 1 / 2 Ar c \~TrZ) [r-J ■ m 

Equations (|B2(I and i|B4() demonstrate that A = (2 In A) 1 / 2 is indeed a very weak function of a when «;(T + )Eo(a) ^> 1, 
thus justifying the constancy of the ratio H\/h s h in the inner optically thick parts of the disk. 



CALCULATION OF THE RADIAL STRUCTURE. 

Numerical solutions for the disk structure based on approximation J7J are obtained using iterative procedure. On 
the first step we adopt some more or less arbitrary initial radial profile for Hi(a), which yields a(a) and T s /j(a) through 
equations © and 113( 1 . This allows us to integrate equations © and (|15fl using temperature profile J7J) and the ideal 
gas law, resulting in distributions of T*(a, z), p(a, z), and P(a, z). With this we can determine new profile of -Hi (a) as 
the vertical height where = 1 at every a. Then we repeat iteration until the convergence is obtained. 

This procedure is generally quite unstable because of the presence of derivative of H\/a with respect to a needed 
for calculation of a. To abate this problem we average each value of Hi(a) over the 5 neighbouring grid points. Also, 
we put a "limiter" in the determination of a, i.e. we limit the relative change of a with respect to its value in the 
previous iteration to be no more that 10%. These simple steps allow us to eliminate spurious oscillations arising in 
the determination of a. We found this procedure to more robust and accurate than the one proposed by Chiang et 
al. (2001) as it allows us to significantly increase radial resolution and obtain very accurate results for the radial disk 
structure without giving rise to parasitic instabilities. 
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